legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
axis.line = element_line(color = "black"),
axis.title.x=element_text(vjust=-0.3, size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times")) +
annotate("text", x = 5, y = 90, label = intensity.r2, parse = T,family="Times", size=3)
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
View(d)
source('~/Desktop/service/figures/service_plot.R')
View(d)
source('~/Desktop/service/figures/service_plot.R')
source('~/Desktop/service/figures/sector_energy_plot.R')
# Packages
################################################
library(doSNOW)
library(foreach)
library(truncnorm)
library(snow)
library(data.table)
library(magrittr)
library(tcltk)
# Variables
conf = 0.05
n.iteration = 10^4.2
lower = F
# data
setwd("/home/blair/Desktop/service/data/world bank")
d = fread("data_format.csv")
d = d[country != "Lesotho"]
# Packages
################################################
library(doSNOW)
library(foreach)
library(truncnorm)
library(snow)
library(data.table)
library(magrittr)
library(tcltk)
# Variables
conf = 0.05
n.iteration = 10^4.2
lower = F
# data
setwd("/home/blair/Desktop/service/data/world bank")
d = fread("data_format.csv")
d = d[country != "Lesotho"]
keep = data.table(d$EG.USE.PCAP.KG.OE, d$EG.EGY.PRIM.PP.KD) %>% na.omit()
x = keep$V1*4.187*10^7/10^9
y = keep$V2
n = length(x)
# Packages
################################################
library(doSNOW)
library(foreach)
library(truncnorm)
library(snow)
library(data.table)
library(magrittr)
library(tcltk)
# Variables
conf = 0.05
n.iteration = 10^4.2
lower = F
# data
setwd("/home/blair/Desktop/service/data/world bank")
d = fread("data_format.csv")
keep = data.table(d$EG.USE.PCAP.KG.OE, d$EG.EGY.PRIM.PP.KD) %>% na.omit()
x = keep$V1*4.187*10^7/10^9
y = keep$V2
n = length(x)
plot(x, y, log = "xy")
plot(y, x, log = "xy")
keep = data.table(d$EG.EGY.PRIM.PP.KD, d$SL.SRV.EMPL.ZS) %>% na.omit()
x = keep$V1*4.187*10^7/10^9
y = keep$V2
n = length(x)
plot(x, y, log = "xy")
r =   nls(log(y) ~ log(100*plnorm(x, meanlog = mu, sdlog = sigma, lower.tail = lower, log.p = FALSE)), start=list(mu = 2, sigma = 2) )
lower = T
r =   nls(log(y) ~ log(100*plnorm(x, meanlog = mu, sdlog = sigma, lower.tail = lower, log.p = FALSE)), start=list(mu = 2, sigma = 2) )
p = predict(r)
beta.est = coef(r)
lower = F
r =   nls(log(y) ~ log(100*plnorm(x, meanlog = mu, sdlog = sigma, lower.tail = lower, log.p = FALSE)), start=list(mu = 2, sigma = 2) )
p = predict(r)
r =   nls(log(y) ~ log(100*plnorm(x, meanlog = mu, sdlog = sigma, lower.tail = lower, log.p = FALSE)), start=list(mu = 10, sigma = 10) )
plot(x, y, log = "x")
r =   nls(y ~ log(100*plnorm(x, meanlog = mu, sdlog = sigma, lower.tail = lower, log.p = FALSE)), start=list(mu = 10, sigma = 10) )
error.func = function(mu, sigma){
predict = log(100*plnorm(x, meanlog = mu, sdlog = sigma, lower.tail = lower, log.p = FALSE))
error = sum((predict - x)^2)
return(error)
}
error.func(2,2)
solution <- nlminb(start = c(1,1),  error.func, lower = c(-Inf, 0), upper = c(Inf, Inf), control = list(iter.max = 100))
solution <- nlminb(start = c(mu = 1,sigma = 1),  error.func, lower = c(-Inf, 0), upper = c(Inf, Inf), control = list(iter.max = 100))
?nlminb
solution <- nlminb(start = c(1, 1),  error.func,
lower = c(-Inf, 0), upper = c(Inf, Inf), control = list(iter.max = 100))
solution <- nlminb(start = c(1, 1),  error.func)
predict = log(100*plnorm(x, meanlog = par[1], sdlog = par[2], lower.tail = lower, log.p = FALSE))
error = sum((predict - x)^2)
error.func = function(par){
predict = log(100*plnorm(x, meanlog = par[1], sdlog = par[2], lower.tail = lower, log.p = FALSE))
error = sum((predict - x)^2)
return(error)
}
error.func(2,2)
solution <- nlminb(start = c(1, 1),  error.func)
solution <- nlminb(start = c(1, 1),  error.func)
beta.est = solution$par
plot(x, y, log = "xy")
points(x.new, p1, col = "red")
points(x.new, p2, col = "red")
p = 100*plnorm(x.new, meanlog = beta.est[1], sdlog = beta.est[2], lower.tail = lower, log.p = FALSE)
x.new = exp(seq(log(1), log(2000), length.out = 200))
p = 100*plnorm(x.new, meanlog = beta.est[1], sdlog = beta.est[2], lower.tail = lower, log.p = FALSE)
p1 = 100*plnorm(x.new, meanlog = lwr$mu2, sdlog = lwr$sigma2, lower.tail = lower, log.p = FALSE)
p2 = 100*plnorm(x.new, meanlog = upr$mu1, sdlog = upr$sigma1, lower.tail = lower, log.p = FALSE)
plot(x, y, log = "xy")
points(x.new, p1, col = "red")
points(x.new, p2, col = "red")
points(x.new, p, col = "red")
plot(x.new, p)
plot(x.new, p, log = "xy")
plot(x, y, log = "xy")
error.func = function(par){
predict = log(100*plnorm(x, meanlog = par[1], sdlog = par[2], lower.tail = lower, log.p = FALSE))
error = sum((predict - log(x))^2)
return(error)
}
solution <- nlminb(start = c(1, 1),  error.func)
solution <- nlminb(start = c(1, 1),  error.func, lower = c(-Inf, 0), upper = c(Inf, Inf), control = list(iter.max = 100))
error = sum((predict - log(y))^2)
return(error)
error.func = function(par){
predict = log(100*plnorm(x, meanlog = par[1], sdlog = par[2], lower.tail = lower, log.p = FALSE))
error = sum((predict - log(y))^2)
return(error)
}
solution <- nlminb(start = c(1, 1),  error.func, lower = c(-Inf, 0), upper = c(Inf, Inf), control = list(iter.max = 100))
beta.est = solution$par
p =  100*plnorm(x, meanlog = beta.est[1], sdlog = beta.est[2], lower.tail = lower, log.p = FALSE)
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
install.packages("rlecuyer")
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
plot(x, y, log = "xy")
points(x.new, p1, col = "red")
points(x.new, p2, col = "red")
points(x.new, p, col = "red")
x.new = exp(seq(log(1e-5), log(10), length.out = 200))
p = 100*plnorm(x.new, meanlog = beta.est[1], sdlog = beta.est[2], lower.tail = lower, log.p = FALSE)
p1 = 100*plnorm(x.new, meanlog = lwr$mu2, sdlog = lwr$sigma2, lower.tail = lower, log.p = FALSE)
p2 = 100*plnorm(x.new, meanlog = upr$mu1, sdlog = upr$sigma1, lower.tail = lower, log.p = FALSE)
plot(x, y, log = "xy")
points(x.new, p1, col = "red")
points(x.new, p2, col = "red")
points(x.new, p, col = "red")
?plnorm
plot(qlnorm(1:100))
plot(qlnorm(seq(0,1, 0.1)))
plot(qlnorm(seq(0,1, 0,001)))
plot(qlnorm(seq(0,1, 0.001)))
error.func = function(par){
predict = qlnorm(y/100, par[1], par[2])
error = sum((predict - log(x))^2)
return(error)
}
solution <- nlminb(start = c(1, 1),  error.func, lower = c(-Inf, 0), upper = c(Inf, Inf), control = list(iter.max = 100))
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
plot(x, y, log = "xy")
points(x.new, p1, col = "red")
points(x.new, p2, col = "red")
points(x.new, p, col = "red")
plot(x.new, p1, col = "red")
points(x.new, p2, col = "red")
points(x.new, p, col = "red")
p = 100*plnorm(x.new, meanlog = beta.est[1], sdlog = beta.est[2], lower.tail = lower, log.p = FALSE)
plot(qlnorm(seq(0,1, 0.001)))
plot(qlnorm(seq(0,1, 0.001)), log = "xy")
error.func = function(par){
predict = qlnorm(y/100, par[1], par[2])
error = sum((predict - x)^2)
return(error)
}
solution <- nlminb(start = c(1, 1),  error.func, lower = c(-Inf, 0), upper = c(Inf, Inf), control = list(iter.max = 100))
beta.est = solution$par
p =  100*plnorm(x, meanlog = 10, sdlog = 1, lower.tail = lower, log.p = FALSE)
plot(x, y, log = "xy")
lines(x, p, col = "red")
p =  100*plnorm(x, meanlog = 1, sdlog = 1, lower.tail = lower, log.p = FALSE)
lines(x, p, col = "red")
p =  100*plnorm(x, meanlog = -5, sdlog = 1, lower.tail = lower, log.p = FALSE)
points(x, p, col = "red")
p =  100*plnorm(x, meanlog = 0, sdlog = 0.2, lower.tail = lower, log.p = FALSE)
plot(x, y, log = "xy")
points(x, p, col = "red")
p =  100*plnorm(x, meanlog = -1, sdlog = 0.2, lower.tail = lower, log.p = FALSE)
points(x, p, col = "red")
p =  100*plnorm(x, meanlog = -1, sdlog = 0.3, lower.tail = lower, log.p = FALSE)
points(x, p, col = "red")
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
source('~/Desktop/service/data/world bank/regression_GDP_intensity.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/sector_intensity_plot.R')
wd = gsub("figures", "bea", dir)
setwd(wd)
wd
wd = gsub("figures", "/data/bea", dir)
setwd(wd)
d = fread("sector_energy_intensity.csv")
test = melt(d, value.name = c("industry", "service"))
test = melt(d, value.name = "industry")
View(test)
View(test)
d = fread("sector_energy_intensity.csv")
d = melt(d, value.name = "industry", value)
d = melt(d)
View(d)
d = melt(d, value.name = "year")
d = fread("sector_energy_intensity.csv")
d = melt(d, value.name = "year")
d = fread("sector_energy_intensity.csv")
d = melt(d, id = year)
d = melt(d, id.vars = "year")
View(test)
d = fread("sector_energy_intensity.csv")
testd = data.table(c(d$industry, d$service),
rep(c("Industry", "Service"), each = length(d$year)
)
text.size = 10
energy_intensity = ggplot(data = d) +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
axis.line = element_line(color = "black"),
axis.title.x=element_text(vjust=-0.3, size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
# export
############################################################
setwd(dir)
gA = ggplotGrob(energy_intensity)
pdf("sector_intensity_plot.pdf", width = 4, height = 3)
grid.arrange( gA )
dev.off()
testd = data.table(c(d$industry, d$service),
rep(c("Industry", "Service"), each = length(d$year))
)
long = data.table(c(d$industry, d$service),
rep(c("Industry", "Service"), each = length(d$year))
)
long = data.table(intensity = c(d$industry, d$service),
sector = rep(c("Industry", "Service"), each = length(d$year))
)
energy_intensity = ggplot(data = long) +
geom_boxplot(aes(x = sector, y = intensity)) +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
axis.line = element_line(color = "black"),
axis.title.x=element_text(vjust=-0.3, size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
source('~/Desktop/service/figures/sector_intensity_plot.R')
energy_intensity = ggplot(data = long) +
geom_boxplot(aes(x = sector, y = intensity)) +
scale_y_continuous("MJ per $", breaks = seq(0, 20,1)) +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
axis.line = element_line(color = "black"),
axis.title.x=element_text(vjust=-0.3, size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
source('~/Desktop/service/figures/sector_intensity_plot.R')
source('~/Desktop/service/figures/sector_intensity_plot.R')
source('~/Desktop/service/figures/sector_intensity_plot.R')
source('~/Desktop/service/figures/sector_intensity_plot.R')
source('~/Desktop/service/figures/sector_intensity_plot.R')
source('~/Desktop/service/figures/sector_intensity_plot.R')
source('~/Desktop/service/figures/sector_intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
setwd(dir)
dir
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
install.packages("egg")
library(egg)
png("intensity_plot.png", width = 6, height = 3,  units = 'in', res = 600)
ggarrange(US, global, widths = 1:2)
dev.off()
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
gtable(gA, gB, widths = c(1,2))
gtable(gA, gB)
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
install.packages("cowplot")
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/intensity_plot.R')
source('~/Desktop/service/figures/sector_energy_plot.R')
dir = here()
wd = gsub("figures", "US_historical", dir)
setwd(wd)
wd = gsub("figures", "data/US_historical", dir)
setwd(wd)
d = fread("sector_energy_intensity.csv")
wd = gsub("figures", "data/US_historical", dir)
setwd(wd)
d = fread("splice_results.csv")
d = d[,c(1, 7:9)]
reshape(d, timevar = year)
reshape(d, timevar = year, "long")
reshape(d, timevar = year, direction = "long")
test = melt(d, id.vars=c("agr_frac", "ind_frac", "srv_frac"))
View(test)
test = melt(d, id.vars=c("year")
text.size = 10
us_sector_time = ggplot(data = long) +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
axis.line = element_line(color = "black"),
axis.title.x = element_blank(),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
# export
############################################################
setwd(dir)
gA = ggplotGrob(us_sector_time)
pdf("us_sectors.pdf", width = 4, height = 3)
grid.arrange( gA )
dev.off()
test = melt(d, id.vars=c("year"))
View(test)
long = melt(d, id.vars=c("year"))
f = function(x){
quantile(x, probs = (c(0.025, 0.5, 0.975)))
}
test = long[, f(value)]
test = long[, f(value), by = variable]
View(test)
test = long[, f(value), by = .(variable, year)]
View(test)
test = long[, as.list(f(value)), by = .(variable, year)]
View(test)
View(test)
result = long[, as.list(f(value)), by = .(variable, year)]
names(result) = c("sector", "year", "lwr", "mid", "upr")
us_sector_time = ggplot() +
geom_ribbon(data = result, aes(x = year, ymin = lwr, ymax = upr, fill = sector)) +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
axis.line = element_line(color = "black"),
axis.title.x = element_blank(),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
source('~/Desktop/service/figures/us_sector_plot.R')
dir = here()
library(data.table)
library(ggplot2)
library(here)
library(gridExtra)
library(grid)
library(gtable)
dir = here()
wd = gsub("figures", "data/eia", dir)
setwd(wd)
dir = here()
wd = gsub("figures", "data/eia", dir)
wd = gsub("figures", "data/eia", dir)
setwd(wd)
e = fread("energy_pc.csv")
conversion = (10^16*1055.06)
e$e_pc*conversion
e$e_pc = e$e_pc*conversion
conversion = (10^16*1055.06)/1000/10^9
e$e_pc = e$e_pc*conversion
library(data.table)
library(ggplot2)
library(here)
library(gridExtra)
library(grid)
library(gtable)
dir = here()
wd = gsub("figures", "data/world bank", dir)
setwd(wd)
d = fread("data_format.csv")
d = d[country != "Lesotho"]
d$EG.USE.PCAP.KG.OE = d$EG.USE.PCAP.KG.OE*4.187/100
wd = gsub("figures", "data/eia", dir)
setwd(wd)
e = fread("energy_pc.csv")
# energy in quad btu
conversion = (10^16*1055.06)/1000/10^9
e$e_pc = e$e_pc*conversion
conversion = (10^12*1055.06)/1000/10^9
library(data.table)
library(ggplot2)
library(here)
library(gridExtra)
library(grid)
library(gtable)
dir = here()
wd = gsub("figures", "data/world bank", dir)
setwd(wd)
d = fread("data_format.csv")
d = d[country != "Lesotho"]
d$EG.USE.PCAP.KG.OE = d$EG.USE.PCAP.KG.OE*4.187/100
wd = gsub("figures", "data/eia", dir)
setwd(wd)
e = fread("energy_pc.csv")
# energy in trillion btu
conversion = (10^12*1055.06)/1000/10^9
e$e_pc = e$e_pc*conversion
plot(e)
